home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 22 / Cream of the Crop 22.iso / os2 / splot170.zip / demo / fftfilt.spt < prev    next >
Text File  |  1996-10-10  |  3KB  |  148 lines
  1.  
  2. #include <splot.h>
  3.  
  4. double *dat;
  5. double *a;
  6. double yint,slope;
  7. int i,j,k,n,n2;
  8. char tmp[80];
  9.  
  10. // datafile name
  11. char *datname = "learn.dat";
  12.  
  13. // filter parameters in hertz
  14. int cutoff = 75;
  15. int rolloff = 25;
  16. int notch = -20;
  17. int width = 10;
  18.  
  19. // data scaling factors
  20. double lpos_scale = 0.234;  // to um
  21. double time_scale = 0.408;  // to ms
  22.  
  23.  
  24. main()
  25.    {
  26.    readdata(datname,dat);
  27.    j = arraydim(dat,0);
  28.  
  29.    // make n next largest power of 2
  30.    n = 1;
  31.    while (n < j) n <<= 1;
  32.    n2 = n * n;
  33.  
  34.    // temp arrays
  35.    double freq[n][3];
  36.    double filt[n][2];
  37.    double lpos[n][2];
  38.    double pow[n/2 +1][2];
  39.  
  40.    // copy lens positions and scale
  41.    for (i = 0; i < j;i++)
  42.       {
  43.       lpos[i][0] = time_scale * i;
  44.       lpos[i][1] = dat[i][3] * lpos_scale;
  45.       }
  46.    // pad end of position array
  47.    for (i = j; i < n;i++)
  48.       {
  49.       lpos[i][0] = time_scale * i;
  50.       lpos[i][1] = dat[j-1][3] * lpos_scale;
  51.       }
  52.  
  53.    // subtract trend
  54.    set(XRANGE,80,280);
  55.    fitline(lpos,0,1,&yint,&slope);
  56.    set(XRANGE,-1e6,1e6);
  57.    for (i = 0; i < n;i++)
  58.       {
  59.       lpos[i][1] = lpos[i][1] - (yint + i * slope);
  60.       freq[i][0] = (double) i / n / time_scale * 1000;
  61.       freq[i][1] = lpos[i][1];   // real part
  62.       freq[i][2] = 0;            // complex part
  63.       }
  64.  
  65.    // transform
  66.    fft(freq,1,2);
  67.    //plotdata(freq,0,1);
  68.    //fft(freq,1,2,-1);
  69.  
  70.    /* build filter */
  71.    int icut,iroll,inotch,iwidth;
  72.    /* convert filt params to integer */
  73.    icut = (int) (cutoff * n * time_scale / 1000);
  74.    iroll = (int) (rolloff * n * time_scale / 1000);
  75.    inotch = (int) (notch * n * time_scale / 1000);
  76.    iwidth = (int) (width * n * time_scale / 1000);
  77.    for (i = 1; i < n / 2;i++)
  78.       {
  79.       filt[i][0] = (double) i / n / time_scale * 1000;
  80.       if (i < icut)
  81.          {
  82.          filt[i][1] = 1.0;
  83.          filt[n-i][1] = 1.0;
  84.          }
  85.       else if (i < icut + iroll)
  86.          {
  87.          filt[i][1] = 1.0 - (double) (i - icut) / iroll;
  88.          filt[n-i][1] = 1.0 - (double) (i - icut) / iroll;
  89.          }
  90.       else
  91.          {
  92.          filt[i][1] = 0;
  93.          filt[n-i][1] = 0;
  94.          }
  95.       if (i >= inotch - iwidth && i <= inotch + iwidth)
  96.          {
  97.          filt[i][1] = fabs((double) (i - inotch) / iwidth);
  98.          filt[n-i][1] = fabs((double) (i - inotch) / iwidth);
  99.          }
  100.       }
  101.  
  102.    // filter data
  103.    for (i = 1; i < n;i++)
  104.       {
  105.       freq[i][1] *= filt[i][1];
  106.       freq[i][2] *= filt[i][1];
  107.       }
  108.  
  109.    // transform back
  110.    fft(freq,1,2,-1);
  111.    for (i = 1; i < n;i++)
  112.       {
  113.       freq[i][0] = lpos[i][0];
  114.       }
  115.  
  116.  
  117.    // plot results
  118.    abox();
  119.    ascale(XAXES,0,300);
  120.    set(AXESCLIP,ON);
  121. #if 1
  122.    ascale(YAXES,-300,100);
  123.    plotdata(lpos,0,1,freq,0,1);
  124.    label(BOTTOM,"Time (ms)");
  125.    label(LEFT,"Lens position (!m!m)");
  126.    text( 8.41,11.91,"Extreme Kink");
  127.    sprintf(tmp,"filter cutoff %d Hz", cutoff);
  128.    text(tmp);
  129. #else
  130.    for (i = 1; i < n;i++)
  131.       {
  132.       lpos[i][1] = lpos[i][1] - freq[i][1];
  133.       }
  134.    plotdata(lpos,0,1);
  135.    label(BOTTOM,"Time (ms)");
  136.    label(LEFT,"Delta Lens position (!m!m)");
  137.    text( 8.41,18.59,"Extreme Kink");
  138.    sprintf(tmp,"filter cutoff %d Hz", cutoff);
  139.    text(tmp);
  140. #endif
  141.    label(TOP,"TH55 drum spin 100");
  142.  
  143.    set(FONTMULT,0.5);
  144.    text(0.5,0.5,CURFILENAME);
  145.  
  146.    }
  147.  
  148.